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I. INTRODUCTION 

Physics of one dimensional (ID) Bose gases attracts 
more and more attention due to both very interesting 
phenomena that appear in these systems and a very 
recent experimental progress toward realization of Tonks 
gases 0, IE • In this paper we provide a unified descrip- 
tion of basic shock waves properties in ID Bose gases in- 
teracting with arbitrary strength. Previous works on this 
subject were limited to either weakly interacting gases 
[E IE IE El or the strongly interacting ones 0] . We 
also describe influence of depleted atoms on dynamics of 
a Bose-Einstein condensate (BEC) shocks. Finally, we 
derive a simple expression for speed of propagation of 
arbitrarily shaped density pulses. 

Sec. |n] describes theoretical basics of an approximate 
hydrodynamic approach. Sec. Ifffl discusses shock-wave 
solutions of this approach on specific examples. Sec. II VI 
describes propagation of arbitrarily shaped pulses. Sec. 
Ivl explains experimental methods used for generation of 
shock structures, while Sec. I VII presents quantum correc- 
tions to shock solutions of the Gross-Pitaevskii equation. 
Finally, Sec. IVIII provides a summary of this paper. 



II. THE MODEL 

We consider delta interacting Bose gas in a ID box. Its 
experimental realization is possible due to recent experi- 
mental progress in trapping of small samples of ultracold 
bosonic atoms in a box- like optical trap Q. Another 
exciting experimental setup may come from a paper of 
Gupta et al [ll|, where successful realization of a ring 
shaped magnetic trap filled with a BEC was recently re- 
ported. 

The Hamiltonian of our system in dimensionless quan- 
tities (see Appendix lAl for units) reads: 
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H = -~22q^ +a2_ J 8{x i -x j ), (1) 
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that there is only one parameter that determines system 
static properties: 
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where p is atomic density normalized to N 3> 1. In the 
limit of 7 — * our system is a Bose-Einstein condensate, 
while when 7 — * +00 it is a Tonks gas, whose proper- 
ties are captured by the Fermi-Bose mapping theorem 
(FBMT) O. To see how big is a BEC-Tonks crossover 
one can look at sound velocity and compare it to BEC 
and Tonks predictions. Exact expression for sound ve- 
locity 01 ( see a l so Appendix is 
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where 6(7) is defined by Lieb and Liniger in 15] . The 
v s /p is depicted in Fig. notice how slowly the system 
enters a Tonks regime. A large BEC-Tonks crossover, 
7 ~ (1,50), provides us motivation for studies of shock 
waves outside BEC H US IE 13 and Tonks HE] limits. 

Many-body solutions of (JTJ at zero absolute tempera- 
ture have been analyzed in classic papers 0, 0] . These 
solutions are based on the Bethe ansatz and allow for an- 
alytical extraction of different static system properties. 
Unfortunately, they are too complicated for description 
of system dynamics. 

To simplify the problem a proper hydrodynamical ap- 
proach can be worked out [la . fl7j | and leads to the fol- 
lowing set of equations 
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where a > is the interaction coupling |l2j ] . It turns out 



P(p) = \p 2 (^(7) - 7^ 



(4) 
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and Vi is an external potential acting on atoms. A nice 
property of these equations is that they exactly repro- 
duce sound velocity © , while the problems with their us- 
age for description of shock propagation (not formation) 
come from their derivation valid in the long-wavelength 
limit, where the quantum pressure (QP) term, 



2 y/P ' 



(7) 



is negligible. 

To get more insight into physics described by JSJ and 
validity of a QP term we look at BEC and Tonks limits. 
In the BEC case, p(p) — > ap and Eqs. 1] 11511 can be 
obtained by time-dependent variational principle applied 
to a product wave function 
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FIG. 1: Sound velocity divided by density. In the BEC limit 
(7 — > 0) it is ^7, while in the Tonks limit (7 — > +00) it equals 
7T. Solid line: ©, dots: Tonks value, dashed line: ^y. For 
units see Appendix 1X1 



,x N ,t) = VN(/)(xiit) ■ ■ ■ 4>(x N ,t), (8) 



III. SHOCK WAVE SOLUTIONS 



j dx\4>(x,t)\ 2 = 1. It results in time-dependent Gross- 
Pitaevskii equation l|29|) . Then substitution of = 
v /7jexp(ix) and v = d x \ into lt2^|) gives a BEC version 
of 11151) . In other terms, p in 11151) is a single particle 
density defined in a many-body theory as 



p(x,t) 



dx 2 ■ ■ ■ dx N \ 1 $(x 1 x 2 , ■ ■ ■ ,x N ,t)\' 
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which equals N\4>(x, t)\ 2 in the BEC limit. A quantum 
pressure term in this limit is rigorously derived [l8| ■ 

In the Tonks regime one has p{p) — ► p 2 w 2 /2, which 
was first found by renormalization group approach in 
[19L l20| . and then used in a number of papers, e.g., 
Oil I2H l22j . Now derivation of © , do not involve any 
sort of product simplification JSJ of a wave-function. In- 
deed, a Fermi-Bosc mapping theorem ^| implies, e.g., 
that if typ(xi, . . . , Xn) is a ground state wave- function of 
noninteracting fermions placed in the same potential as 
a Tonks gas, then a ground state Tonks wave-function is 
\^ f{xi, ■ ■ ■ , xn)\- Therefore, even a simple description 
of a Tonks gas based on the single particle density © 
involves N orthogonal single particle orbitals instead of 
a single one, 4>(x), used in the BEC case. 

In Tonks limit the QP term was shown to lead to un- 
physical density oscillations 0, l23| . From knowledge 
that the quantum pressure term is present in the BEC 
limit 7^0 and absent in Tonks regime 7 — > +00 it is 
clear that Eqs. 1 1151) can not describe shock propaga- 
tion for arbitrary 7. Nonetheless, shock formation from 
density perturbations that are initially wide can be suc- 
cessfully done, and will be described below. In this case 
the QP term is unimportant roughly up to shock for- 
mation. Due to lack of theoretical concepts for getting 
exact time-dependent solutions for the system of inter- 
est, we consider future experimental results as the best 
verification of our calculations based on hydrodynamic 
equations. 



The quantum pressure term (jTJ is important only when 
density changes occur on length scales smaller then the 
characteristic length given by 
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In a BEC limit £(p) is called a healing length, and we 
propose to use the same name regardless of 7. 

Now we want to solve Eqs. (1 ll5l(i[) . Assuming that 
the perturbation under consideration is, at least initially, 
broad compared to £, the quantum pressure term in © 
can be neglected. Following standard methods [24| one 
gets 



p{x,t) = f 
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To proceed further with analytical calculations we use 
the following relations: 
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where 7 C ~ 14.5 from the requirement that both expan- 
sions match at j c . Expression i|13fl is extracted from 
Lieb's observation [lj| that for 7 <~ 10 

v / 7 -7 3/ 7(2 7 r). (15) 
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where the first term corresponds to the Fermi length of a 
non-interacting Fermi gas - a result that may be expected 
from the FBMT [13. 

Rewriting implicit solution to the form 

g(x, t) = f (x - Coo Woo (g) t) , 

where Coo = tt/?o is a background sound velocity of a 
Tonks gas, one gets from <|11[) and i|14|l 



2g + sg 2 - 1 
(l + sp) 2 (l + s) 
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Taking (|16|l as an initial density profile, the implicit 
shock-wave solution reads: 



g(x,t) 
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cosh[(x - CooWoc^t) / a] 
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FIG. 2: Dotted line: initial density profile (lib! . Solid line: 
solution of Eqs. 14151 with density at t — given by 1161 and 
velocity field determined from 11 1 21 . Dashed line: implicit so- 
lution (list multiplied by po- Upper plot: subsequent density 
profiles correspond to t — 0, 1.155, 2.02. Lower plot: t = 2.6. 
Time of shock creation equals t s — 2.03, iV = 499, 70 = 30, 
a = 2,1 = 50, 77 = 0.05. For units see Appendix 1X1 



Since Q15fl complicates further calculations, we expanded 
it into a series around yfiy/(2ir) equal to zero. Expression 
(|14fl is taken directly from Lieb and Liniger paper [l5j - 
see Eqs. (3.32) and Appendix lAl 

To test above approximations we compare the shock 
wave solutions that neglect the QP term and use ap- 
proximate expression for J ^§^> to the full numerical 
solution of hydrodynamic equations (|4I5I6[I . Initially we 



determine velocity field from 1|12|) and choose 



p(x,0) = po 



VPo 



COsh(x/cr) 



V > 0, 



(16) 



with a ^> £(p) (calculations for 77 < can be easily re- 
peated). The background density po is found from the 

normalization condition: J, dx p(x, 0) = N, where we 
have assumed that periodic box has boundaries at ±1. 
For well localized perturbations being of interest from 
now on one gets po = 2 (i+ V a) ' Fi nau y 5 it is convenient to 
define relative density 



P 

Q = — ■ 

Po 

Now we consider separately 7 > j c and 7 < j c cases. 



A. 7 > 7c case 



The healing length is found from (|10|l and Eqs. (3.32) 
of subjected to rescalings of Appendix 1X1 



(17) 
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Although, the explicit form of g{x, t) can not be ex- 
tracted analytically, important properties of pulse dy- 
namics can be determined. First of all, impulse ampli- 
tude, I + 77, is constant in time. Second, speed of impulse 
maximum for any 77 > equals 

V(r)) = CooWoo(l + J?), 

and applicability of this expression is not limited to 
l/cosh(x/a) 2 perturbations only - Sec. IIVI Third, the 
width of the impulse at given density is constant during 
propagation- a property that we missed in earlier papers 
[jj • Fourth, impulse tails propagate roughly with the 
background sound velocity equal to c 00 w 00 (l), while the 
impulse maximum moves with the speed CooWoo(l + ??)• 
Since w(g) monotonically increases with g, the impulse 
deforms its shape so that a shock wave front forms, i. e., 
\d x p(x,t)\ = +00 at one point. 

Time and position of shock-wave creation can be ex- 
tracted from equations [24| : 



d p x(p) = , d 2 p x(p)=0 



(19) 



Their solution gives: density p s at which density profile 
becomes locally vertical, and time t s at which this occurs: 



6s = 



3 + 4s?7 + 9s - y/9{l + s) 2 + 16sV 
6s 



ts = 



Vy/V (1 + SQs) 3 



4Coo {Qs - + V - Qs 



(20) 



(21) 



These expressions close to a Tonks gas limit take a simple 
form: 
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Comparison of analytical solution l|18(l and full numer- 
ical one based on hydrodynamic equations (|4I5I6|I is pre- 
sented in Fig. |3 As easily seen, there is a good agreement 
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between hydrodynamical solution and a shock-wave one 
until the moment of shock creation. Then discrepancy 
increases due to appearance of density oscillations trig- 
gered by the QP term (JJJ neglected in derivation of (|18H . 
Since presence of the QP term for large 7, e.g., 7 > j c , is 
questionable, it is an interesting open question whether 
density oscillations in the form presented in Fig. [21 sur- 
vive or not for any 7 > 7 C considered here. 



B. 



7 < 7c case 



The healing length reads 



V^ap 
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(22) 



derived by extraction of p(p) from (|15() and subsequent 
Taylor expansion of HlOjl . Naturally, the first term in 122(1 
corresponds to the BEC healing length. 

The implicit shock-wave solution has the form: 

g(x,t) = f (x - cowo(g)t) , (23) 

where cq — ^/apo is the background speed of sound in 
the limit of 7 — > 0. Combining and l|13l) one gets 



Wq(q) = 3y/g-2 
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(24) 



To verify accuracy of l|23|l we have simulated dynamics for 
the system with 70 = a/po = 5, i. e., for 70 outside the 
BEC mean-field regime (see inset of Fig. QJ, but small 
enough to stay clearly within 7 < 7 C parameter range. 
Agreement between full numerical solution of Eqs. 1 1101) 
and implicit shock solution is good before shock forma- 
tion. Close to time of shock creation discrepancies show 
up in the oscillatory region missed in 1(23)) . Qualitatively, 
the plot that presents these results is the same as Fig. [21 
Since dynamics of initial density profile is qualitatively 
the same as in the 7 > 7 C case, we notice that impulse 
maximum moves now for any rj > with velocity 

V(r?) = c w (l + n)- 

The explicit solution of shock equations i|19|) can be 
found in the limit of 70 — * 0, but still arbitrary 77 > 0: 

1 + 7?+ V(l +y)(9 + y) 

Ps = -, , 



t. 



3c (p s - 1)^/1 + 77 - p s 



Interestingly, for a gaussian like initial impulse discussed 
in [5j explicit expressions for p s and t s in this limit were 
beyond the reach. Finally, we note that the QP term for 
7 — ► is certainly present, but it is hard to say whether 
it survives in JJjJ for any 7 < j c . Since the QP term 
affects shock dynamics strongly, an experiment should 
clarify this uncertainty. 



IV. PROPAGATION OF PULSES OF 
ARBITRARY SHAPE 

In this section we describe some general properties of 
shock- wave solutions 1)11112)1 . To this aim we consider 
density profiles satisfying 



p(x,t) = f(x-W(p)t) = f((). 



(25) 



Suppose that there are n extrema (minima, maxima), 
placed on a background density po, in the initial density 
profile p(x,0). Let us denote positions of these extrema 
at t = as Xi(0) and set pi — p(xi(0),0). It means that 
dxf(x)\ Xt (o) = and <9^/( x )\xi(o) 7^ 0- At t > one gets 
from (|2SJl 



d x p{x, t) 



df(Q 



1 



+ dpW{p f_m t 



which obviously equals zero iff £ — Xi{Q). It implies 
that Xi(t) — Xi(0) + W(pi)t. Furthermore, we have 
dxP(x,t)\ Xi{t) = dlp{x, 0)1^(0) + 1 0, which proves that 
extrema do not change into saddle points in the course 
of time evolution. This moans that velocity of impulse ex- 
trema equals W(pi), where pi = p(xi(t),t) = p(xi(0),0). 

Since W(p) is shape independent, we conclude that 
for any 7 density extrema propagate with constant am- 
plitude and speed equal to W(pi). This speed is well 
approximated as 



< 7 < 7 C 



CoW (pi/p ), 



7 > 7c : CooWooipi/po), 



(26) 



(27) 



for 7's a little off j c . 

Finally, it is important to stress that above statements 
concerning the amplitude and velocity of density extrema 
are valid at least as long as the implicit solution works, i. 
e., up to the time of shock creation. After that time, they 
are valid in regions far enough from shock structures. 



V. EXPERIMENTAL REALIZATION 

As discussed in 0, 0, experimental creation of 
matter-wave packets undergoing shock-wave dynamics is 
straightforward. The idea is to cool atoms in an addi- 
tional external potential created by a well-detuned laser 
beam, and then suddenly turn the beam off. The term 
suddenly, means on time scale much smaller than time of 
sound propagation through the perturbation. 

The initial density of atoms resulting from an external 
laser potential Vi is determined as follows. Combining 
density p and velocity field v = d x \ into 4> = v / pexp(ix) 
the Eqs. (gEl) can be rewritten to the standard form 



id t 



1 



dl(f> + Vi{x)(/> + p(p)<l>, 



(28) 
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which time-independent version reads 

--^+y ( (i)^,i(^ = ^ 

with p being a chemical potential. Assuming that the 
laser potential induces density changes on length scales 
larger than the healing length Q1U[I. one determines an 
initial density profile from 

Vi{x)+n{p) = p, 

with p found from J dx p = N . Naturally, velocity field 
equals zero in an initial state considered here. 

A straightforward generalization of results from [H ITo| 
shows that if the initial density profile is p(x,0) = 
po + h(x) with max|ft,(a;)| <C po then p{x,t) — po + 
h (x - v s (po)t) /2 + h (x + v s (p Q )t) /2, with v s (p ) being 
background sound velocity 10) ■ Therefore, there are left 
and right moving perturbations. Due to mirror symme- 
try, we can forget about the left moving one once the 
pulses are well separated. The right-moving pulse is de- 
scribed by po + h(x — v s {po)t)/2 as long as t <C t s , since 
even arbitrarily small density perturbations experience a 
shock deformation after long enough evolution. 

In BEC and Tonks limits qualitatively the same split- 
tin g p rocess takes place even when max|/i(a:)|/yOo ~ 1 
[El ll(|. When the initial density profile is p(x, 0) = 
po + h(x), the right moving pulse is described by p(x, t) w 
Po+f{x— W(p)t) with W(p) being the same as in (|23|) and 
f(x) — h(x)/2. We checked by solving hydrodynamic 
equations (|4I5[) that above described splitting process 
works qualitatively the same way for any 7, as it works 
for 7 — * 0,+oo 0, EJ- This observation explains how 
we imagine practical generation of matter-wave pulses 
described in previous sections. 



VI. APPLICABILITY OF THE 
GROSS-PITAEVSKII MEAN-FIELD APPROACH 
TO SHOCK PROBLEMS 

In the BEC limit system dynamics in the mean-field 
approximation (JHJ is described by the following version 
ofEq. 

id t 4> = -\d 2 x (t> + Vi{x)<t> + aN\^\ 2 4> = H GP 4>, (29) 

where J dx \(f>(x, t)\ 2 — 1. This is the Gross-Pitaevskii 
(GP) equation, which can be rigorously derived by dif- 
ferent means - see, e.g., [T^.l25|. 

To go beyond the mean-field approximation we em- 
ploy a second quantization formalism, i. e., we transform 
Hamiltonian Q with an additional potential term, Vi(x), 
into the form 

H = Jdx¥(-^ + V l (x) + ~m^)*, (30) 



where ^ is a field operator. In the number conserving 
Bogolubov approach |2{| it reads 

* = a (/) + 

with So annihilating a particle from a condensate, and 
annihilating a particle from modes orthogonal to the 
condensate one. Naturally, <fi is a condensate mode, i. 
e., it is an eigenstate of a single particle density matrix, 
(^(x') 1 ^(x)), to the highest eigenvalue (A — N): 

J dx' (&(x')V(x))<f>(x') = \ (f>(x). 

The single particle density of atoms, (W (x) 1 ^ (x)) , equals 
(ala Mx)\ 2 +J2H(x)\ 2 , (31) 

k 

where the first term accounts for condensate density 
while the second one corresponds to density of depleted 
atoms. In further calculations we need both Uk and Vk, e. 
g., for finite temperature calculations. Indeed, density of 
atoms at sufficiently low temperature T (Xo(T)/N ~ 1) 
equals [25l | 

(dla )\^)\ 2 +J2( h l^)(uk\u k ) +J2^ k + 

k k 

where {b\b k ) = [exp(Lu k /(k B T) - To get u k , v k 

and to k one first finds a condensate mode from the time- 
independent GP equation, Hgp4> = M^j an d then con- 
structs the matrix 

( H GP + aNQ\4>\ 2 -p aNQ4> 2 \ 

^ -aNQ*(j)* 2 -H GP - aNQ*\(f>\ 2 + p J ' 

where Q is a projector to space orthogonal to a conden- 
sate mode: Qtp = ip—<j}{<f>\ij;). The modes and frequencies 
are found from the eigen equation 



( u k > 






]=UJm {v k k) 







while their dynamics is captured by 

4( uk )=c( uk ), 

dt\v k J V Vk J 

where p is set to zero in £, and the calculation of opera- 
tor C at every time step requires simultaneous solution of 
the time-dependent Gross-Pitaevskii equation. Expres- 
sions for densities of atoms l|31|l and l|32|l are unchanged 
in a time dependent case. For a list of u, v modes prop- 
erties, and details of the number-conserving Bogolubov 
approach see |25j. 

It is important to realize what are limitations of the 
Bogolubov approach. To simplify notation we discuss 
T = case now (the results for T ^ can be easily ob- 
tained) . The Bogolubov method fails in the following two 
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situations. First, the global breakdown, shows up when 
total number of depleted atoms, i. e., J dxj^k |^fe(^)| 2 , 
becomes comparable to number of atoms in the system. 
This limitation comes directly from derivation of Bogol- 
ubov Hamiltonian |25| . Second, local breakdown, hap- 
pens when density of depleted atoms becomes compara- 
ble to condensate density on a length scale of the order 
of the healing length: 



/ dx'N\<j>(x')\ 2 ~ / cfe'VMx')| 2 , (33) 
Jx-fi Jx-f£ 



with / = 0(1). To clarify (|33|l . we notice that in the 
Bogolubov approach an impact of noncondensed atoms 
on condensed ones is neglected. Indeed, to get Bogol- 
ubov dynamics one solves a time-dependent GP equation 
and then puts this solution into equations of motion of 
noncondensed atoms. In this way noncondensed parti- 
cles are influenced by condensed ones but not other way 
round. Next order corrections take into account change 
of condensate density as a result of repulsive interactions 
between condensed and noncondensed clouds, and lead 
to modifications of the GP equation. The correction to 
the right hand side of 129|) is qualitatively of the form 
Vdc P i</> (Eqs. (95,96) in 25]), where 



Vdepl (x) = 2a \vk{x)Y 



(34) 



so that Vdepi (x) can be regarded as an external potential 
acting on a condensate. Extension of changes of conden- 
sate density induced by an external, localized potential, 
is given by the healing length [l^. Therefore, impact of 
depleted atoms on a condensate is best given by com- 
parison of the new term, Vd op i averaged over a few heal- 
ing lengths, to previously dominating term in the GP 
equation, aN\(f>(x)\ 2 , averaged on the same set of points. 
It gives the condition (|33|l . Appearance of such a local 
breakdown of Bogolubov approach will be shown after 
breaking of a dark shock density profile resulting in soli- 
ton production. 

Finally, let us comment when condensate depletion 
might affect mean-field predictions. Since density of de- 
pleted atoms depends on the product of a and N only, it 
is unaffected by the limit 



TV 







aN = fixed, 



while the condensed part in total atomic density grows 
as (ajao) ~ N 1*311320 . It means that in this limit, be- 
ing to some extent possible with the help of Feshbach 
resonances |26|. density of depleted atoms does not af- 
fect density measurements. The exception might happen 
when there is dynamical instability leading to fast in- 
crease of depleted fraction. Then sooner or later density 
of depleted atoms will start to have an important impact 
on condensate dynamics. 

In a generic experimental situation, however, a frac- 
tion of depleted atoms is non-negligible. To compare a 
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FIG. 3: Density of a condensate, (a,Qao)\(j)(x)\ 2 , at t = 0,0.6: 
(a) and (c). Density of depleted atoms Q^ fe=0 \ v k{ x )\ 2 ) a t 
t = 0,0.6: (b) and (d). Plot (e): comparison between conden- 
sate density (solid line) and density of depleted atoms (dashed 
line) in the right-moving shock structure; condensate density 
is rescaled for presentational purposes. At t = laser poten- 
tial Vi was oc — exp(— 2x' 2 ) and the system was in a ground 
state; Vi(t > 0) = 0. The number of noncondensed atoms 
equals ~ 320 during whole evolution. Total number of atoms 
is 1.5 • 10 5 . Other parameters: I = 15, aN = 7500, T = 0. 
Time of shock creation obtained from 5] is t s = 0.36. For 
units see Appendix 1*51 



mean- field prediction to quantum many-body dynamics, 
we have chosen a particular set of experimentally acces- 
sible parameters - see Appendix 1*51 for details. 

In the following we discuss evolution of shock struc- 
tures resulting from white (rj > 0) and dark (77 < 0) 
initial density profiles, where r\ stands for relative per- 
turbation amplitude: as in Eq. 1*16*1 . Since corrections 
to mean field equations are very different in these two 
cases we discuss them separately. 

White shock waves: We have done different cal- 
culations changing not only the size of an initial per- 
turbation but also aN, N, etc. In all cases we observe 
that total depletion of a condensate negligibly increases 
from beginning of time evolution to appearance of well 
developed shock structures. Density of depleted atoms 
around shocks becomes at most 0(1) larger than deple- 
tion density far away from perturbations. Therefore, if 
density and distribution of depleted atoms was such that 
the Bogolubov approach was initially applicable, it will 
also work during shock creation and propagation. 

The results for experimentally relevant choice of pa- 
rameters are depicted in Fig. [3J As easily seen from 
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differences in scaies between Figs. Eta), (c) and Figs. 
Etb), (d), the corrections to the mean- field result are mi- 
nor. Repulsion between condensed and depleted atoms, 
results in localization of depleted atoms around conden- 
sate density minima: Fig. OJe). This was predicted by 
us in 0], however, the amount of depleted atoms, turns 
out to be insufficient to fill condensate ripples. 

We have also done finite temperature calculations ac- 
cording to i|32|) . At T 7^ not only quantum depletion 
but also thermal one shows up, which leads to increase 
of total number of depleted atoms in comparison to the 
T = situation. We have considered low temperatures, 
i. e., T < 30nK (see Appendix |Ej| , to see whether ther- 
mal effects qualitatively affect shock dynamics. The only 
difference with respect to T = case is that total deple- 
tion is larger by a factor of O(10) for T ~ 30nK, which 
is still not enough for getting significant corrections to 
BEC shocks from depleted atoms. 

These calculations suggest that the Gross-Pitaevskii 
equation captures correctly physics of white shock waves 
in a BEC. Finite temperature effects does not seem to 
destroy the qualitative picture of shock formation and 
propagation given by the mean-field approach. 

Dark shock waves: Dark density profiles undergo 
two distinct stages of evolution. In the first one, shock 
waves form, while in the second one shock density profile 
breaks into a train of solitons that move, according to the 
mean-field approach, with constant velocity and without 
changes of shape. The same was independently observed 
numerically in [27j and theoretically in [2^ . Theoretical 
predictions in both papers are based solely on the basis 
of mean-field equations. 

As in white shocks calculations, we have fixed the pa- 
rameters to those experimentally relevant- see Appendix 
IbI The formation of shock wave structures is depicted 
in Figs. Eta), (b), (c), (d). Now a rear impulse edge 
sclf-steepens instead of a front one. This comes from 
the fact that now impulse tail moves faster than impulse 
center, i. e., density minimum. In the case considered 
in Fig. EJ a shock profile breaks before complete sepa- 
ration. Obviously, if the impulse would be more shallow 
the breakdown would happen later. As easily seen, the 
correction to total density coming from depleted atoms 
can be neglected at this stage of time evolution. 

After shock breakdown, solitons are produced as de- 
picted in Fig. Etc)- Then density of depleted atoms in- 
creases by orders of magnitude. It results in significant 
corrections to total density of atoms and causes local, i. 
e., around soliton minima, breakdown of the Bogolubov 
approach according to l|33ll . Notice that still fraction of 
depleted atoms, less than 0.7% in Fig. E[f), is so small 
that the global condition of Bogolubov approach appli- 
cability is well satisfied. As depicted in Fig. [5]both peak 
density of depleted atoms at soliton minimum (dots) and 
total depletion of a condensate (squares) follow approx- 
imately power law increase from the moment of soliton 
train formation. 

It is easy to predict qualitatively what are correc- 
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FIG. 4: Density of a condensate, (a^ao) \(f>(x)\ 2 , at t = 0,0.06: 

(a) and (c). Density of depleted atoms Q^fctfo l Wfc (- E )| 2 ) a ^ * ~ 
0,0.06: (b) and (d). Plot (e): black thick line - total density 
of atoms, dashed line - condensate density; both curves are for 
right-moving dark-shock-originated soliton train at t = 0.75. 
Plot (f): density of depleted atoms at t = 0.75. External laser 
potential at t = was oc exp(— 2a; 2 ) and the system was in 
a ground state; Vi(t > 0) = 0. The number of noncondensed 
atoms equals ~ 320 at t = and ~ 1000 at t = 0.75. The 
total number of atoms equals 1.5 ■ 10 5 . Other parameters: 
I = 15, aN = 7500, T = 0. Time of shock creation obtained 
from [f| is t s = 0.086. For units see Appendix iFll 



tions to a Bogolubov result coming from depleted atoms. 
The next order corrections to the Bogolubov approach 
come from inclusion of repulsive interactions between 
condensed and non-condensed atoms. It is qualita- 
tively accounted for by introducing a new potential term, 
Eq. (|34p. into the Gross-Pitaevskii equation. It means 
that soliton structures in condensate density will become 
wider, so that depleted atoms will have more space to dis- 
tribute themselves inside solitons. Therefore, increase of 
peak density of depleted atoms will certainly slow down. 
This observation is supported by nonperturbative calcu- 
lation of Dziarmaga 29] , who found that total density of 
atoms at soliton minimum approaches background den- 
sity for large times. 

Such a fast increase in density of depleted atoms 
around soliton minima is rather a generic phenomenon. 
Indeed, it was first theoretically discussed for phase- 
imprinted solitons in harmonic traps |30l l3l| as a mech- 
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FIG. 5: Introducing notation: p - peak density of depleted 
atoms in the deepest soliton moving to the right (Fig. 
dN - total number of depleted atoms. Dots: lnp(£); dotted 
line: a linear fit to circles, with slope given by 2.29. Squares: 
\n[dN(t) — dN(0)]; dashed line: linear fit with slope equal to 
2.55. The data is plotted from t = 0.24, i.e., a moment of 
the deepest soliton creation. For the parameters: see Fig. 2] 
while for units see Appendix ITU 



anism for fast disappearance of solitons in the Hannover 
experiment |32| . The experimental time scale for soli- 
ton disappearance was found both theoretically l3~H 
and experimentally to be O(10ms). Our calcula- 
tion, based on the same value of nonlinear parameter 
aN = 7500, gives the same estimation for instant when 
corrections to dynamics of dark-shock-originated solitons 
become important. Indeed, the time instant 0.75 in Fig. 
Efe), corresponds to ~ 12ms- see Appendix iBl for a unit 
of time. 

To end this section we recall that predictions presented 
here apply to a ID BEC system. There was recently 
an attempt to experimentally verify applicability of the 
mean- field approach to shock phenomena |33l | . The con- 
clusion was that the ripples in a white shock wave front 
are not filled with depleted atoms. Our present work sup- 
ports this observation assuming that ID results are qual- 
itatively correct in 2D rapidly rotating array of vortices 
on top of which a white shock wave front propagates. 



VII. SUMMARY 

This paper presents a complete hydrodynamical de- 
scription of basic properties of shock-waves in a delta 
interacting ID Bose gas. Our predictions can be directly 
verified experimentally. In principle, a proper experimen- 
tal setup should answer at what Lieb-Liniger parameter 
7 the quantum pressure term starts to lead to unphysi- 
cal results. Notice that density oscillations in front of a 
shock wave are present due to this term for small enough 
7, while for large 7 they are absent in exact solution and 
shocks propagate in a very different way |Toj . In fact, 
a change of this kind in shock dynamics can be a nice 
experimental signature of a BEC-Tonks crossover. 



We also discussed quantum many-body corrections to 
Gross-Pitaevskii shock-wave solutions. This way we clar- 
ified the role of depleted atoms in shock dynamics. This 
important point was missed in previous studies of BEC 
shock waves jfj, |6J, |7J, |8j, |9( . 

From experimental point of view, tools for verification 
of our predictions seem to be available either right now or 
in the nearest future. Indeed, a first experiment on sound 
propagation in a BEC was done years ago Q ■ Since we 
do not discuss time-of-flight measurements the best com- 
parison between our theory and experiment should rely 
on in situ measurements |34| . The most prospective ex- 
perimental setup for observation of shock dynamics in 
different 7 parameter ranges is probably provided by a 
box-like optical trap Q , where values of 7 < 1 approach- 
ing the BEC-Tonks crossover (Fig. ^) , have been already 
achieved. Other interesting experimental setups include 
atom chips |35| and circular atom waveguides The 
latter one, being best suited for studies of shock colli- 
sions. 

From theoretical side, there are at least two possible 
interesting extensions of this work. First of all, one can 
analyze shock dynamics in quasi-lD trapping geometries, 
where the system is a 3D waveguide with adjustable har- 
monic transverse confinement. This can be easily done 
using results from [3(|. Indeed, instead of Eq. (JSJ one 
can consider Eq. (5) of 36] and repeat subsequent cal- 
culations. In the limit of tight transversal confinement, 
i. e., for a ID Bose gas, both expressions for n(p) lead to 
physically identical results. The second interesting exten- 
sion of this work includes studies of possible outcomes of 
a single shot density measurement. In fact, it should be 
stressed that predictions based on hydrodynamical and 
Bogolubov approaches apply to averages over different 
experimental measurements done on the system prepared 
many times in the same quantum state. Needless to say , 
averages may differ from a single shot outcomes - see [371 ] 
for an illustrative example. 

I would like to acknowledge discussions with Zbyszek 
Karkuszewski. This work was started in the Institute for 
Theoretical Physics in Hannover. I'm grateful to both the 
Alexander von Humboldt Foundation for support of this 
work in Germany, and to the US Department of Energy 
for support of this research at the LANL. 



APPENDIX A: UNITS IN A BOX 

Since this paper uses extensively results of 0, 0] it 
is useful to link dimensionless quantities used by us to 
those of 

The eigen equation of a ID Hamiltonian expressed in 
terms of dimensional quantities denoted by primes is 

1,2 E^+^'E 5 ^-^)*-^ ( A1 ) 



2m 



to get dimensionless Hamiltonian Q one introduces: 
x'i = Xilo, E' = Eh 2 /(mlQ), A' = ah 2 /(mlo), where Iq 
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is an arbitrary length scale. Consideration of a time- 
dependent Schrodinger equation leads to t' = tnil^/h. 

The dimensionless Hamiltonian used by Lieb and Lin- 
iger has the form 

d 2 



i<3 



and can be obtained from I|A1I) after rescalings: x\ — Xilo, 
E' = Eh 2 /{2mll), A' = ch 2 /{ml ), t' = t2ml 2 /h. 

Therefore, there are the following relations between 
our dimensionless quantities and Lieb and Liniger ones 
marked by LL: 



7 = 111 



/' : 



(J>LL 



VLL 

2 : 



where v is sound velocity. Finally we note, that there is 
a misprint in the first line of expression (1.4) providing 
sound velocity 0|. There should be —2^/de/d'y (as used 
in Sec. [H| instead of —jde/d'y. The second line of (1.4) 
in [Ijl is correct. 



where J dx \4>(x)\ 2 = 1 and l± is a variational parameter. 
Substituting IB 211 into (|B1|I one gets 



£[<j>] = J dx 
which supports the following ID GP equation 



l\\ 2 J_ 
2 + 2lY 



(B3) 



id^ = + y + aN\<P\ 2 ^ ( B-l ) 



with aN = g^^, / (2tt1 2 _) . To determine l±, we substitute 
the Thomas-Fermi solution of 1|B4(I . 



aN\<j)(x)\ 



2 _ 1 /3aN\ 2/3 x 2 



(B5) 



into (|B3|I getting £(l±), which minimization leads to 



APPENDIX B: A BOX APPROXIMATION OF A 
3D SYSTEM IN THE BEC LIMIT 

We aim at getting a ID box approximation of a full 
3D BEC gas confined in a 3D harmonic potential. The 
system in the mean-field approximation satisfies the fol- 
lowing 3D Gross-Pitaevskii equation 



iti 



a*' 



V' 2 *' + (mtu 2 x' 2 /2 + muj 2 , r' 2 /2W 

2m 



where c<j_l 3> oj, g' 3D = AnTi 2 a sc N/m with a sc being inter- 
atomic scattering length, and J d 3 x' \ty'\ 2 — 1. First, we 
rescale all quantities using the harmonic oscillator units 
in the x' direction. Denoting by symbols without primes 
dimensionless quantities, 



ui 



(x', y>, z<) = (x, »,Ji,*' = « (^) 3/4 
V tyiuj V n / 

we arrive at the following 3D GP equation 

= + {x2/2 + AV / 2 )* + 93d\*\ 2 *, 

where A is trap aspect ratio ujj_/lo while g^p — 
Aira sc N \J muj/fi. This equation can be derived from the 
following energy functional: 



d A x 



1 '— - 12 , f„2/o , \2„2/oM,T,|2 , 93D |,j t |4 



|V*| 2 + (2;72 + AV/2)|1f 



(Bl) 

To proceed further we assume a simple form of a varia- 
tional wave-function, 



^(x, r) = 4>(x) exp 



1 



2l\) 0Fk' 



(B2) 



(^k) 2/3 = (2507T 2 /9) 1 /3 (A 2 /i-l). 

In the limit of g^o — > one gets l±_ = 1/VA, which is a 
noninteracting cloud width. In a typical experiment, a 
cloud width is much larger than 1/VX H|, so 
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V 250tt 2 



A 3 / 5 ' 



For a cigar shaped cloud one can assume, e.g., u> — 2n ■ 
10Hz, wi = 2vr • 569Hz. Taking also N = 1.5 • 10 5 , a sc = 
5.2-10~ 9 m (scattering length of 87 Rb in the \F = 2, m F = 
2) state), m = 89.91 x 1.66 • 10~ 27 kg (atomic mass of 
87 Rb), one can find that the units of length and time are 
~ 3.4/um and ~ 16ms, respectively. For these parameters 
aN w 7500 and the unit of temperature, Tiuj/kB, equals 
0.48nK. By using all these results it is easy to transform 
dimensionless plots from Sec. I VII into dimensional ones. 

A box approximation of dimensionally reduced har- 
monically trapped cloud described above is the following. 
We place exactly the same number of atoms in the box as 
in the quasi-lD configuration described by (|B5|I . We as- 
sume that, in the absence of external laser potential, den- 
sity in a box extending from [— I, I] is exactly the same as 
at a center of a harmonic trap: [3/(4v / 2)] 2 ^ 3 /(aA r ) 1 / 3 = 
1/(20, which gi ves 1 = (aA r ) 1/3 (2/3) 2 / 3 . 

For parameters defined above 2-1 w 30 (about 0.1mm), 
which can be compared to the Thomas-Fermi size of the 
harmonically trapped cloud equal here to ~ 44 (about 
0.15mm). The box-like approximation leads to results be- 
ing in good qualitative agreement with calculations done 
in a harmonically trapped case. It concerns both length 
and time scales of shock dynamics. 
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